In [31]:
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import pypandoc
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from sklearn import linear_model
import matplotlib.pyplot as plt
from matplotlib import rcParams
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.nonparametric.smoothers_lowess import lowess
import numpy as np
import scipy.stats as stats
from sklearn.preprocessing import PolynomialFeatures
Exercise 9¶
9 a) Produce scatterplot matrix for all variables¶
In [2]:
auto = pd.read_csv('Auto.csv')
In [3]:
fig = px.scatter_matrix(auto)
In [4]:
fig.show()
9 b) compute matrix of correlations¶
In [5]:
correlation_cols = set(auto.columns)
correlation_cols.remove('name')
In [6]:
corr = auto[correlation_cols].corr()
In [7]:
corr.style.background_gradient(cmap='coolwarm')
Out[7]:
9 c) multiple linear regression with mpg as response¶
In [9]:
auto['horsepower'] = pd.to_numeric(auto['horsepower'],errors='coerce')
In [10]:
auto = auto.dropna(axis=0,how='any')
In [11]:
input_cols = correlation_cols
input_cols.remove('mpg')
In [12]:
auto.shape
Out[12]:
In [13]:
input_cols
Out[13]:
In [14]:
X = auto[input_cols]
In [15]:
y = auto[['mpg']]
In [16]:
X = sm.add_constant(X)
In [17]:
model = sm.OLS(y,X).fit()
In [18]:
predictions = model.predict(X)
In [19]:
model.summary()
Out[19]:
9 c i) Is there a relationship between predictors and response?¶
Yes, there is a relationship between the predictors and the response based on the r-squared value of 0.821, which shows that 82% of the mpg response can be explained by the predictors in the model. Some predictors don't have as much statistical significance to the response though, based on their p-values.
9 c ii) Which predictors appear to have statistically significant relationship to response¶
Displacement, Weight, Origin and Year appear to have a statistically significant relationship to the response based on their p-values.
9 c iii) What does the coefficient for year suggest?¶
If all other predictors held constant, the mpg increases by 0.75 for every increase in year.
9 d) Diagnostic plots of linear regression fit¶
In [21]:
residuals = model.resid
fitted = model.fittedvalues
smoothed = lowess(residuals,fitted)
Residual Plot¶
In [22]:
fig, ax = plt.subplots()
ax.scatter(fitted, residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Residuals')
ax.set_xlabel('Fitted Values')
ax.set_title('Residuals vs. Fitted')
ax.plot([min(fitted),max(fitted)],[0,0],color = 'k',linestyle = ':', alpha = .3)
Out[22]:
This plot shows a non-linear relationship between the response and the predictors
Normal Q-Q Plot¶
In [23]:
sorted_student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
sorted_student_residuals.index = model.resid.index
sorted_student_residuals = sorted_student_residuals.sort_values(ascending = True)
In [24]:
df = pd.DataFrame(sorted_student_residuals)
df.columns = ['sorted_student_residuals']
df['theoretical_quantiles'] = stats.probplot(df['sorted_student_residuals'], dist = 'norm', fit = False)[0]
rankings = abs(df['sorted_student_residuals']).sort_values(ascending = False)
fig, ax = plt.subplots()
x = df['theoretical_quantiles']
y = df['sorted_student_residuals']
ax.scatter(x,y, edgecolor = 'k',facecolor = 'none')
ax.set_title('Normal Q-Q')
ax.set_ylabel('Standardized Residuals')
ax.set_xlabel('Theoretical Quantiles')
ax.plot([np.min([x,y]),np.max([x,y])],[np.min([x,y]),np.max([x,y])], color = 'r', ls = '--')
Out[24]:
This plot shows that the data shows a mostly normal distribution but slightly right skewed
Scale-Location Plot¶
In [25]:
student_residuals = model.get_influence().resid_studentized_internal
sqrt_student_residuals = pd.Series(np.sqrt(np.abs(student_residuals)))
sqrt_student_residuals.index = model.resid.index
smoothed = lowess(sqrt_student_residuals,fitted)
In [26]:
fig, ax = plt.subplots()
ax.scatter(fitted, sqrt_student_residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('$\sqrt{|Studentized \ Residuals|}$')
ax.set_xlabel('Fitted Values')
ax.set_title('Scale-Location')
ax.set_ylim(0,max(sqrt_student_residuals)+0.1)
Out[26]:
This plot shows that the residuals are randomly spread across the range of the predictors.
Residuals vs. Leverage Plot¶
In [27]:
student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
student_residuals.index = model.resid.index
df = pd.DataFrame(student_residuals)
df.columns = ['student_residuals']
df['leverage'] = model.get_influence().hat_matrix_diag
smoothed = lowess(df['student_residuals'],df['leverage'])
sorted_student_residuals = abs(df['student_residuals']).sort_values(ascending = False)
In [28]:
fig, ax = plt.subplots()
x = df['leverage']
y = df['student_residuals']
xpos = max(x)+max(x)*0.01
ax.scatter(x, y, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Studentized Residuals')
ax.set_xlabel('Leverage')
ax.set_title('Residuals vs. Leverage')
ax.set_ylim(min(y)-min(y)*0.15,max(y)+max(y)*0.15)
ax.set_xlim(-0.01,max(x)+max(x)*0.05)
plt.tight_layout()
cooksx = np.linspace(min(x), xpos, 50)
p = len(model.params)
poscooks1y = np.sqrt((p*(1-cooksx))/cooksx)
poscooks05y = np.sqrt(0.5*(p*(1-cooksx))/cooksx)
negcooks1y = -np.sqrt((p*(1-cooksx))/cooksx)
negcooks05y = -np.sqrt(0.5*(p*(1-cooksx))/cooksx)
ax.plot(cooksx,poscooks1y,label = "Cook's Distance", ls = ':', color = 'r')
ax.plot(cooksx,poscooks05y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks1y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks05y, ls = ':', color = 'r')
ax.plot([0,0],ax.get_ylim(), ls=":", alpha = .3, color = 'k')
ax.plot(ax.get_xlim(), [0,0], ls=":", alpha = .3, color = 'k')
ax.annotate('1.0', xy = (xpos, poscooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, poscooks05y[-1]), color = 'r')
ax.annotate('1.0', xy = (xpos, negcooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, negcooks05y[-1]), color = 'r')
ax.legend()
plt.show()
This plot shows a potential leverage at leverage > 0.175 point which is well outside the other points leverage value
Do the residual plots suggest any unusually large outliers?¶
The residual plots do not show many unusually large outliers, but the leverage plot does identify one observation with large leverage
9 e) Fit linear regression models with interaction effects¶
In [29]:
X = auto[input_cols]
In [32]:
x_interaction = PolynomialFeatures(2, interaction_only=True, include_bias=False).fit_transform(X)
In [33]:
interaction_df = pd.DataFrame(x_interaction, index=auto.index, columns = ['cylinders','displacement','horsepower','weight','acceleration','year','origin',
'cylinders:displacement','cylinders:horsepower','cylinders:weight','cylinders:acceleration',
'cylinders:year','cylinders:origin','displacement:horsepower','displacement:weight',
'displacement:acceleration','displacement:year','displacement:origin','horsepower:weight',
'horsepower:acceleration','horsepower:year','horsepower:origin','weight:acceleration',
'weight:year','weight:origin','acceleration:year','acceleration:origin','year:origin'])
In [34]:
y = auto[['mpg']]
In [35]:
interaction_df = sm.add_constant(interaction_df)
In [36]:
model = sm.OLS(y,interaction_df).fit()
In [37]:
predictions = model.predict(interaction_df)
In [38]:
model.summary()
Out[38]:
In [39]:
model.pvalues[model.pvalues < 0.05]
Out[39]:
In the interaction model above with all of the interactions included the displacement:weight, displacement:origin, and acceleration:origin interactions all appear to be significant. This shows that although origin does not appear to be significant on its own it is important in relationship to displacement and acceleration so it should be included in the model
9 f) Try a few different transformations of the variables¶
In [40]:
auto.columns
Out[40]:
In [43]:
auto['sqrt_displacement']=auto['displacement']**(1/2)
auto['log_displacement']=np.log(auto['displacement'])
auto['displacement_squared']=auto['displacement']**(2)
In [44]:
fig = px.scatter(auto, x='sqrt_displacement', y='mpg')
fig.show()
In [45]:
fig = px.scatter(auto, x='log_displacement', y='mpg')
fig.show()
In [46]:
fig = px.scatter(auto, x='displacement_squared', y='mpg')
fig.show()
In [47]:
fig = px.scatter(auto, x='displacement', y='mpg')
fig.show()
From the plots above for displacement, sqrt(displacment), log(displacment) and displacement^2, it appears log and square root show more linear relationship with mpg than displacement and dispalcement^2.
Exercise 10¶
10 a) Fit multiple regression model to predict sales using price, urban and us¶
In [57]:
carseats = pd.read_csv('Carseats.csv')
In [58]:
carseats['Urban'] = carseats.Urban.map(dict(Yes=1, No=0))
carseats['US'] = carseats.US.map(dict(Yes=1, No=0))
In [59]:
carseat_input_cols = ['Price', 'Urban', 'US']
In [60]:
X = carseats[carseat_input_cols]
In [61]:
y = carseats['Sales']
In [62]:
X = sm.add_constant(X)
In [63]:
model = sm.OLS(y,X).fit()
In [64]:
predictions = model.predict(X)
In [65]:
model.summary()
Out[65]:
10 b) Interpretations of each coefficient in the model¶
When the other two predictors are held constant, when price is increased by $1, sales dicreases by 0.0545 sales. The binary (yes/no) urban predictor is not statistically significant to sales. Stores in the US sell 1200 more units than non US stores, with all other predictors held constant
10 c) Write model in equation form:¶
Sales = 13.0435 - 0.545 Price - 0.0219 Urban + 1.2 * US + $\epsilon$ where Urban and US are binary (Yes=1 and No=0)
10 d) For which predictors can you reject null hypothesis¶
We can reject null hyptohesis for Price and US, as their P-values are less than 0.05
10 e) Fit smaller model based on response to previous question¶
In [66]:
carseat_input_cols = ['Price', 'US']
In [67]:
X = carseats[carseat_input_cols]
In [68]:
y = carseats['Sales']
In [69]:
X = sm.add_constant(X)
In [70]:
model = sm.OLS(y,X).fit()
In [71]:
predictions = model.predict(X)
In [72]:
model.summary()
Out[72]:
10 f) How well do the models in (a) and (e) fit the data?¶
Both have R-squared values of 0.239, meaning they only explain 24% of the response in the prediction.
10 g) Using model from (e), obtain 95% confidence interval for the coefficients¶
From model summary in (e): price: [-0.065, -0.044]; US: [0.692, 1.708]
10 h) Is there evidence of outliers or high leverage observations in model from (e)?¶
In [73]:
residuals = model.resid
fitted = model.fittedvalues
smoothed = lowess(residuals,fitted)
Residual Plot¶
In [74]:
fig, ax = plt.subplots()
ax.scatter(fitted, residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Residuals')
ax.set_xlabel('Fitted Values')
ax.set_title('Residuals vs. Fitted')
ax.plot([min(fitted),max(fitted)],[0,0],color = 'k',linestyle = ':', alpha = .3)
Out[74]:
Normal Q-Q Plot¶
In [75]:
sorted_student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
sorted_student_residuals.index = model.resid.index
sorted_student_residuals = sorted_student_residuals.sort_values(ascending = True)
df = pd.DataFrame(sorted_student_residuals)
df.columns = ['sorted_student_residuals']
df['theoretical_quantiles'] = stats.probplot(df['sorted_student_residuals'], dist = 'norm', fit = False)[0]
rankings = abs(df['sorted_student_residuals']).sort_values(ascending = False)
fig, ax = plt.subplots()
x = df['theoretical_quantiles']
y = df['sorted_student_residuals']
ax.scatter(x,y, edgecolor = 'k',facecolor = 'none')
ax.set_title('Normal Q-Q')
ax.set_ylabel('Standardized Residuals')
ax.set_xlabel('Theoretical Quantiles')
ax.plot([np.min([x,y]),np.max([x,y])],[np.min([x,y]),np.max([x,y])], color = 'r', ls = '--')
Out[75]:
Scale-Location Plot¶
In [76]:
student_residuals = model.get_influence().resid_studentized_internal
sqrt_student_residuals = pd.Series(np.sqrt(np.abs(student_residuals)))
sqrt_student_residuals.index = model.resid.index
smoothed = lowess(sqrt_student_residuals,fitted)
fig, ax = plt.subplots()
ax.scatter(fitted, sqrt_student_residuals, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('$\sqrt{|Studentized \ Residuals|}$')
ax.set_xlabel('Fitted Values')
ax.set_title('Scale-Location')
ax.set_ylim(0,max(sqrt_student_residuals)+0.1)
Out[76]:
Residuals vs. Leverage Plot¶
In [77]:
student_residuals = pd.Series(model.get_influence().resid_studentized_internal)
student_residuals.index = model.resid.index
df = pd.DataFrame(student_residuals)
df.columns = ['student_residuals']
df['leverage'] = model.get_influence().hat_matrix_diag
smoothed = lowess(df['student_residuals'],df['leverage'])
sorted_student_residuals = abs(df['student_residuals']).sort_values(ascending = False)
In [78]:
fig, ax = plt.subplots()
x = df['leverage']
y = df['student_residuals']
xpos = max(x)+max(x)*0.01
ax.scatter(x, y, edgecolors = 'k', facecolors = 'none')
ax.plot(smoothed[:,0],smoothed[:,1],color = 'r')
ax.set_ylabel('Studentized Residuals')
ax.set_xlabel('Leverage')
ax.set_title('Residuals vs. Leverage')
ax.set_ylim(min(y)-min(y)*0.15,max(y)+max(y)*0.15)
ax.set_xlim(-0.01,max(x)+max(x)*0.05)
plt.tight_layout()
cooksx = np.linspace(min(x), xpos, 50)
p = len(model.params)
poscooks1y = np.sqrt((p*(1-cooksx))/cooksx)
poscooks05y = np.sqrt(0.5*(p*(1-cooksx))/cooksx)
negcooks1y = -np.sqrt((p*(1-cooksx))/cooksx)
negcooks05y = -np.sqrt(0.5*(p*(1-cooksx))/cooksx)
ax.plot(cooksx,poscooks1y,label = "Cook's Distance", ls = ':', color = 'r')
ax.plot(cooksx,poscooks05y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks1y, ls = ':', color = 'r')
ax.plot(cooksx,negcooks05y, ls = ':', color = 'r')
ax.plot([0,0],ax.get_ylim(), ls=":", alpha = .3, color = 'k')
ax.plot(ax.get_xlim(), [0,0], ls=":", alpha = .3, color = 'k')
ax.annotate('1.0', xy = (xpos, poscooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, poscooks05y[-1]), color = 'r')
ax.annotate('1.0', xy = (xpos, negcooks1y[-1]), color = 'r')
ax.annotate('0.5', xy = (xpos, negcooks05y[-1]), color = 'r')
ax.legend()
plt.show()
From the Normal Q-Q plot there don't appear to be any outliers and from the Residuals vs. Leverage plot there don't appear to be any leverage points